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Abstract 

Homology models of mammalian voltage-gated sodium (Nay) channels based on the crystal structures of the bacterial 
counterparts are needed to interpret the functional data on sodium channels and understand how they operate. Such 
models would also be invaluable in structure-based design of therapeutics for diseases involving sodium channels such as 
chronic pain and heart diseases. Here we construct a homology model for the pore domain of the Navl.4 channel and use 
the functional data for the binding of //-conotoxin GIIIA to Navl.4 to validate the model. The initial poses for the Navl.4- 
GIIIA complex are obtained using the HADDOCK protein docking program, which are then refined in molecular dynamics 
simulations. The binding mode for the final complex is shown to be in broad agreement with the available mutagenesis 
data. The standard binding free energy, determined from the potential of mean force calculations, is also in good 
agreement with the experimental value. Because the pore domains of Nayl channels are highly homologous, the model 
constructed for Navl.4 will provide an excellent template for other Nayl channels. 
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Introduction 

Voltage-gated sodium (Nay) channels are responsible for 
initiation and propagation of action potential, which is essential 
for the activity of excitable cells such as neurons, heart and muscle 
cells [1]. Nay channels are involved in many physiological 
activities throughout the body, which also makes them potential 
drug targets for related disorders such as cardiac and neuropathic 
diseases [2-5] . Because of their prominence, much effort has gone 
into understanding the structure and function of Nay channels. 
For a very long time, no molecular structures were available for 
Nay channels, and the low resolution structures obtained from 
cryo-electron microscopy [6] were not very informative. In the 
absence of crystal structures, indirect methods such as studying the 
effect of mutations on hgand binding [7-12] have been the main 
source of information on Nay channels. 

This situation has changed dramatically with the recent 
determination of the crystal structures of several bacterial Nay 
channels [13-17]. As in potassium channels, where solution of the 
bacterial K""" channel KcsA [18] started an intense period of 
examination of the structure-function relations, we expect a similar 
thing to happen in Nay channels. Already, there have been many 
computational studies of the bacterial Nay channels, investigating 
their ion permeation [19-27] and gating [28,29] properties, as well 
as ligand binding [30,31]. There are, as yet, no crystal structures 
for the mammalian Nay channels. Unlike in potassium channels, 
the tetrameric symmetry is lost when going from the bacterial to 
the mammalian Nay channels. Therefore, solution of the 
structures of the mammalian Nay channels is likely to be more 



drflicult and may not follow soon. This leaves construction of 
homology models from the bacterial ones as the most viable route 
for making progress. Again this is not as straightforward as in 
potassium channels because, besides the loss of the tetrameric 
symmetry, the critical selectivity filter is not conserved either 
between the bacterial and the mammalian Nay channels. 
Nevertheless, the bacterial Nay channels provide a reasonable 
scaffold for construction of homology models for the mammalian 
ones, and such models can be constrained and validated using the 
large amount of functional data that have been accumulated over 
several decades [1]. Initial attempts in this regard include a 
docking study of tetrodotoxin and anesthetics binding to the pore 
domain of the Nay 1.4 channel [32], and an MD study of Na"'"/K"'' 
selectivity in Nayl channels [33]. 

Nay channels are targets for many toxins, which bind with high 
afifmity to various sites on the channel protein, disabling its normal 
function. For example, tetrodotoxin, saxitoxin, and //-conotoxins 
bind to the channel vestibule and block the pore [1,34]. This has 
been exploited in many experimental studies, where the known 
toxin structures were used to probe the pore domain of Nay 
channels [35,36]. Because ;U-conotoxins are the only peptide toxins 
that block Nay channels, there has been a lot of interest in their 
properties. The first /i-conotoxin to be characterised was jl- 
conotoxin GIIIA (/i-GIIIA), which selectively binds to the Nay 1.4 
channel [37,38]. Its structure was determined from NMR by 
several groups [39-41]. A great deal of mutagenesis and functional 
studies have been performed on the Navl.4-/z-GIIIA complex to 
determine the binding mode and identify the key residues involved 
in binding [42-55]. For this reason. Nay 1.4— ^-GIIIA complex 
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provides a unique system for testing the liomology models of 
mammalian Nay channels that are constructed from the bacterial 
ones. 

Here, we construct a homology model for the pore domain of 
the Nay 1.4 channel based on the crystal structure of NayAb. The 
stability of the Nay 1.4 model is checked via molecular dynamics 
(MD) simulations. We then create a model for the Navl.4— ^- 
GIIIA complex using the docking program HADDOCK [56,57], 
followed by refinement in MD simulations. The proposed binding 
mode for the Nay 1.4— //-GIIIA complex is shown to give a 
satisfactory account of the available mutagenesis data. As a 
dynamical test of the complex model, we have determined the 
binding free energy of //-GIIIA from its potential of mean force 
(PMF), which is again found to be in good agreement with the 
experimental value. The same computational methodology has 
previously been used to study toxin binding to potassium channels, 
where its ability to yield accurate protein-ligand complexes and 
binding free energies has been demonstrated [58-64]. The present 
study extends it to sodium channels, where much more work 
remains to be done. 

Methods 

Modeling of Nav1.4 Channel 

We construct a homology model for the pore domain of Nay 1.4 
using the crystal structure of NayAb (PDB ID: 3RVY) [13]. The 
sequence of Nayl.4 is taken from the Uniprot database (P15390). 
As shown in Figure 1, the sequences for the four domains of 
Nayl.4 are all quite different from that of NayAb. This makes 
homology modeling of Nay 1.4 not so straightforward as has been 
noted earlier [32]. Here, we have aligned the critical (DEKA) 
residues in the selectivity filter with the corresponding (EEEE) 
residues in 3RVY. A gap is placed between the E and W residues 
in domain II of the selectivity to align the conserved W residues. 
Then, we do multiple ahgnments between S5-pore-S6 sequence 
for all four domains of Nayl.4 and NayAb using the program 
ClustalW [66]. The final ahgnments obtained for the P1-SF-P2 
sequences are shown in Figure 1, which are the same as those in 
[15] but differ from the ones in [32] in the placement of the gap in 
domains I, III, and IV, which is after the W residues in the 
selectivity filter. 

While acceptable alignments are obtained for the pore domain, 
it is much harder to do the same for the S5-P1 and P2-S6 linker 
sequences in the turret because there are no good templates, and 
much less data are available on their function. Therefore, we have 
neglected these linker sequences in our current model of Nayl.4. 



While the S5-P1 linker faces the pore, it does not appear to be 
involved in binding of /i-GIIIA, hence our results are unlikely to be 
affected by its absence. A 3D model of the channel is created using 
Modeller [67] by threading the aligned Nayl.4 sequence for each 
domain on a corresponding domain of 3RVY. In order to refine 
the model and check its stability, we have performed MD 
simulations of the Nayl.4 model in a membrane environment. For 
this purpose, we have used the protocols established in previous 
MD simulations of ion channels [68,69]. The Nay 1.4 model is 
embedded in a lipid bUayer consisting of 153 POPE molecules in 
the x-y plane and solvated with a 0.1 M NaCl solution. Extra 
counter ions are included in the system to neutralize it where 
necessary. The system is then equilibrated in MD simulations in 
several stages. First the protein is frxed and the system is 
equilibrated with pressure coupling until the correct water and 
lipid densities are obtained. At that point, the x and y dimensions 
of the simulation box are fixed and pressure coupling is applied 
only in the z direction (the box size is about 95 x96 x84 A'^). In the 
second stage, the restraints on the protein atoms are relaxed 
gradually by first reducing those on the side chain atoms from 
A: = 30 kcal/mol/A^ to 0 in 3 ns. Finally, the backbone atoms are 
relaxed in a similar manner. The resulting system is run for 0. 1 /Js 
to check the stability of the model. The rmsd of the backbone 
atoms of Nay 1 .4 formed a plateau after the first few ns which 
remained stable throughout the 0. 1 /Js of MD simulations with an 
average value of 0.6 A, confirming its stability. 

Modeling of Nav1.4-/i-GIIIA Complex 

The conotoxin /(-GIIIA is a 22-residue peptide with the 
sequence RDCCTOOKKC KDRQCKOQRC CA with an 
amidated Ala at the C-terminal. It has three Arg, four Lys and 
two Asp residues with a total charge of +6«. The NMR structure of 
//-GIIIA is shown in Figure 2A (PDB ID: ITCJ) [41]. The side 
chains of the four basic residues involved in binding of /i-GIIIA to 
Nayl.4 are indicated explicitiy in the figure. The structure of //- 
GIIIA is stabilized by three disulfide bridges between C3-C15, 
C4-C20, and C10-C21. However, as seen from the superposition 
of the NMR snapshots in Figure 2B, the peptide has quite a bit of 
flexibility around the hinges at the C3 and CIO residues. For this 
reason, the residues 1 1-22 forming the binding interface will be 
used in the rmsd calculations when we check for distortion of the 
toxin during umbrella sampling MD simulations. 

The initial poses for the Nayl.4-/;-GIIIA complex are 
generated by docking the //-GIIIA structure to the Nayl.4 model 
using the program HADDOCK [56,57]. In previous studies of 
toxin binding to potassium channels, HADDOCK has given very 
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Figure 1. Alignment used in homology modeling of rNa^ 
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1.4. The DEKA residues in the four domains forming the selectivity filter (SF) are 
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Figure 2. NWIR structures of //-GIIIA. (A) Structure of fi-G\\\A with the pore blocl<ing residues Kl 1, R13 and K16 pointing downward. Three disulfide 
bridges and the C3-K10 hydrogen bond stabilizing the structure are indicated explicitly. (B) Superposition of the ten NMR structures demonstrating 
the flexibility of the N-terminal residues 1-9 around the C3 and CIO hinges. 
doi:1 0.1 371 /journal.pone.01 05300.g002 



good results, reducing the time needed for refinement witli MD 
[59-63,70]. In order to get an adequate sampling of the side chain 
orientations, we use all ten NMR conformers of ^-GIILA in 
ensemble docking. Because there are no well-known binding 
motifs for Nayl channel blockers — like the pore inserting Lys in 
potassium channel blockers — ^we have considered several possibil- 
ities for restraints in HADDOCK. To facilitate comparisons with 
the mutation data and simplify interpretation of the results, we use 
a single restraint in each docking study. The EEDD and DEKA 
ring of residues are the potential sites on the channel for using 
restraints. However, the mutation data indicates that the EEDD 
residues play a much more important role in binding of /i-GIIIA 
than the DEKA residues [47] . Therefore, only the EEDD ring is 
used as a restraint site in the following docking studies. The 
potential restraint sites on the toxin include the residues Rl, Kl 1, 
R13, K16, and R19, which are identified in mutagenesis 
experiments (see Table 1). Separate docking studies are performed 
for each of these residues and the EEDD ring as a restraint. One 
thousand conformations are generated from each docking, and the 
top hundred selected on the basis of scoring are analyzed via 
clustering to find the dominant binding mode. In choosing among 
the five sets of complexes generated, we have also used the fact 
that /(-GIIIA blocks the channel. This criterion is satisfied only for 
the complexes in the set with the R 1 3-EEDD restraint, which also 
have the highest scores. The complex structure in the top-ten of 
this set that has the most toxin-channel contacts is selected for 
refinement in MD simulations (see Table 2). 

The complex structure obtained from docking is aligned with 
the channel model in the membrane, and the coordinates of the 
toxin are transferred to the channel model. The protocol used in 
equilibrating the channel protein is also used for the complex, with 



the channel protein and toxin being relaxed simultaneously. The 
equilibrated system is run for another 50 ns to check its stability 
and generate trajectory data for analysis. During this MD 
simulation, a small restraint with A: = 0.1 kcal/mol/A" is applied 
to channel backbone atoms to preserve the integrity of the channel 
but no restraints are imposed on the toxin. The system is found to 
be well-equilibrated from the start, and therefore, the trajectory 
data are used for the analysis of the complex structure. 

MD Simulations and PMF Calculations 

All MD simulations are performed using version 2.7 of NAMD 
[71] with the CHARMM36 force field [72]. An NpT ensemble is 
used with periodic boundary conditions. Pressure is kept at 1 atm 
and temperature at 300 K using Langevin coupling with damping 
coefficients of 5 ps~ . Lennard-Jones interactions are switched off 
smoothly within distance of 10-13.5 A. Electrostatic interactions 
are computed without truncation using the particle-mesh Ewald 
algorithm. A time step of 2 fs is employed in all MD simulations. 
The trajectory data is saved at 1 ps intervals but the reaction 
coordinate is written at every time step in umbrella sampling 
simulations. 

The PMF for dissociation of //-GIIIA from Nay 1.4 is 
constructed using umbrella sampling MD simulations. As the 
method has been described in detail previously [58,63], we give 
only a brief description here. The reaction coordinate is chosen as 
the distance between the center of mass (COM) of the channel 
protein and the COM of the toxin along the channel axis. Initially 
28 umbrella windows are created along the channel axis at 0.5 A 
intervals using steered MD with a force constant k = 30 kcal/ mol/ 
A^ and pulling speed v = 5 A/ ns. Six more windows are added 
subsequently to ensure that the toxin has reached the bulk region 
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Table 1 


. Effect of the mutations in ^u-GIIIA on the IC50 values for binding 


to Na„1.4. 






Res. 


Mut. 


IC55„(wt) 


IC5o(mut) 


Ref. 






(nM) 


IC5„(wt) 




Rl 


A 


175 


5.0 


[42] 


R1 


A 


11 


154 


[48] 


Rl 


A 


30.7 


17.35 


[53] 


Rl 


K 


175 


2-6 


[42] 


Rl 


Q 


26.6 


1.8 


[43] 


D2 


A 


175 


0.50 


[42] 


D2 


N 


26.6 


0.53 


[43] 


06 


P 


26.6 


0.79 


[43] 


07 


P 


26.6 


0.79 


[43] 


K8 


A 


175 


2-6 


[42] 


K8 


A 


11.2 


22 


[55] 


K8 


Q 


26.6 


7.2 


[43] 


K9 


A 


11.2 


0.36 


[55] 


K9 


Q 


26.6 


1.65 


[43] 


Kll 


A 


175 


2-6 


[42] 


Kll 


A 


30.7 


10.4 


[53] 


Kll 


A 


11.2 


29 


[55] 


Kll 


Q 


26.6 


0.77 


[43] 


D12 


A 


175 


0.28 


[42] 


D12 


A 


11 


0.65 


[48] 


D12 


N 


26.6 


0.67 


[43] 


R13 


A 


175 


229 


[42] 


R13 


A 


11.7 


210 


[46] 


R13 


A 


30.7 


642 


[53] 


R13 


K 


175 


2-6 


[42] 


R13 


K 


11.7 


16.8 


[46] 


R13 


Q 


11.7 


121 


[46] 


R13 


D 


11.7 


764 


[46] 


014 


D 


11 


31.9 


[48] 


K16 


A 


175 


2-6 


[42] 


K16 


A 


11 


24.3 


[48] 


K16 


A 


30.7 


81.4 


[53] 


K16 


Q 


26.6 


2.23 


[43] 


017 


P 


26.6 


0.97 


[43] 


017 


P 


11 


14.3 


[48] 


018 


K 


26.6 


0.53 


[43] 


R19 


A 


175 


25 


[42] 


R19 


A 


30.7 


199 


[53] 


R19 


A 


11.2 


577 


[55] 


R19 


K 


175 


1.0 


[42] 


R19 


Q 


26.6 


2.68 


[43] 


The residue and its mutation are listed in the first two columns. The experimental IC50 values for the wild type are given in the third column and the ratio IC5o(mut)/ 
ICsolwt} in the fourth column. The data are collected using different experimental systems {oocytes, mammalian cells or lipid bilayers} and varying ionic strengths, which 
partly explains the variation In the measured IC50 values. 
doi:10.1371/journal.pone.0105300.t001 



signalled by a flat PMF. The same force constant of /: = 30 kcal/ avoid numerical instabilities in construction of the PMF. We have 

mol/A^ is used in umbrella sampling simulations, which is found included two more windows between the windows 4-5 and 22-23, 

to be optimal for toxins of this size [58]. The overlaps of where this condition is not satisfied. The reaction coordinates 

distributions between the neighboring windows should be > 5% to collected from the simulations are unbiased and combined using 
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the weighted histogram analysis method [73]. Umbrella sampling 
simulations are continued until the convergence of the PMF is 
assured from block data analysis of the PMF data. 

The binding constant is determined by integrating the PMF, 
W(z), along the z axis [58,63] 



-W{z)lk,T 



where the Z\ and Zj are the initial and final points in the PMF, and 
nF? is the average cross sectional area of the binding pocket as 
explored by the COM of the toxin, which is determined from its 
transverse fluctuations. The value of R is obtained from restraint- 
free MD simulations of the Navl.4-/i-GIIIA complex as 0.79 A. 
Finally, the standard binding free energy of //-GIIIA is obtained 
from the binding constant using 



Gy,= -kBT\n{K,^CQ) 
where Co is the standard concentration of 1 M. 



(2) 



Results and Discussion 

Binding Mode of the Nav1.4-^-GIIIA Complex 

Snapshots of the Navl.4-^-GIIIA complex obtained from 
docking and MD simulations are shown in Figure 3 (A PDB file 
giving the coordinates of the complex structure is provided in File 
S2). The five basic residues on the toxin that make contacts with 
the acidic residues on the channel are indicated explicitly. To 
provide a more quantitative picture of the binding mode, we have 
calculated the average N-O distances between the interacting 
residues from the 50 ns MD trajectory data (Table 2). The 
distances obtained from HADDOCK are also shown for 
comparison. Of the nine contact pairs listed, HADDOCK has 
got six of them within 0.3 A of the MD result and one of them 
within 1.5 A. Only in two cases, the discrepancy is larger than 5 A. 
Considering that the binding mode is rather complex involving 
multiple partners for some of the residues, this is an excellent result 
from HADDOCK. It also explains why the complex structure 
obtained from HADDOCK has equilibrated so quickly in MD 
simulations. 

In order to validate the Navl.4-^-GIIIA complex, we compare 
the binding mode characterized in Figure 3 and Table 2 with the 
mutagenesis data summarized in Table 1. The toxin residue R13 



makes the most contacts with the channel residues and hence is 
expected to play a major role in binding, which is in good 
agreement with experiments (Table 1). The R13A mutation 
causes more than 200-fold increase in the IC50 ratio. From Eq. 
2, this corresponds to more than 3 kcal/mol change in the binding 
free energy. In previous studies of binding of ShK toxin to Kyi 
channels, [59-61], we have shown that neutralizing a charged 
toxin residue in contact with channel residues costs about 2 kcal/ 
mol in binding free energy, provided the binding mode is 
preserved after the mutation. However, if the binding mode of 
the mutated toxin is substantially different from that of the wild 
type, there could be larger changes in the binding free energy. We 
have, therefore, performed docking calculations for jl- 
GIIIA[R13A] and found that the binding mode has completely 
changed with Rl inserting in the Nay 1.4 pore. This highlights the 
importance of checking the binding mode when interpreting 
unusual results in alanine scanning experiments. Other mutations 
of R13 also reduce the affinity, including the conservative 
mutation R 1 3K, presumably because a Lys residue cannot sustain 
the multiple contacts R13 makes. R13 is followed in importance 
by the residues R19, K16, and Kll, each of which makes tight 
contacts with the acidic residues on the channel. For these 
residues, the observed loss of affinity, when they are mutated to 
Ala, are mosdy in the range expected from switching off the charge 
(Table 1). The quality of the charge contact between K8 and 
D 1 248 is not as good as the others because K8 is on the flexible N- 
terminal part of the toxin, which fluctuates around the CIO hinge. 
As a result, the K8-D1248 average distance and its sigma are 
substantially larger compared to the other charge contacts. Again 
this is consistent with the experimental observation that K8A 
mutation causes less aflhnity loss relative to the other contacts listed 
in Table 2 (Table 1). Further evidence for the relative strength of 
the individual interactions will be presented when we discuss the 
persistence of contacts during dissociation of ^u-GIIIA from 
Navl.4. 

The only residue whose mutation to Ala reduces the affmity of 
//-GIIIA but does not appear in the proposed binding mode is Rl. 
As seen from Figure 3, Rl is on the opposite side of the binding 
interface and does not make any contacts with the channel 
residues (the N-O distance for the closest acidic residue on the 
channel is more than 9 A). We have seen in studies of potassium 
channel toxins that mutation of an Arg residue could still result in 
a different binding mode even though it is not directly involved in 
binding [59,61]. To see if a similar thing happens in //-GIIIA, we 
have docked //-GIIIA[R1A] to Navl.4. The resulting binding 



Table 2. List of interacting residues in the Nav1.4-//-GIIIA complex. 



/^GIIIA 


Na„1.4 


Dock. 


MD Aver. 


K8-N^ 


01248-0, 


2.6 


4.0±1.2 


Kll-N^ 


D1 241-0, 


7.8 


2.7±0.4 


K11-N^ 


Dl 532-0, 


2.7 


2.6±0.1 


RI3-N2 


E403-O, 


2.6 


2.7±0.4 


RI3-N1 


E758-02 


3.1 


2.8±0.3 


RI3-N2 


DI532-O2 


9.0 


2.7±0.2 


K16-N^ 


E758-0, 


2.7 


2.7±0.1 


K16-N^ 


DI24I-O2 


2.6 


2.7±0.1 


RI9-N2 


D762-02 


2.6 


2.7±0.4 


The average N- 


-0 distances obtained from HADDOCK and MD simulations 


are given in the third and fourth columns {in units of A). 



doi:l 0.1 371/journal.pone.Ol 05300.t002 
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Figure 3. Binding mode of the Nav1.4-/i-GIIIA complex. For clarity, domains I and III (top) and II and IV (bottom) are shown separately. All the 
important interactions between the channel and toxin residues are indicated explicitly. 
doi:1 0.1 371 /journal.pone.01 05300.g003 



mode is found to be substantially different from that of ;U-GIIIA, 
which indicates that the reduction in affinity is likely to be caused 
by a change in the binding mode rather than loss of a charge 
interaction with the channel residues. 

The binding mode emerging from this study provides a 
complete block of the fairly large vestibule of Nay 1.4. This is 
illustrated in Figure 4, which shows a bottom-view of the complex. 
The three basic residues, R13, K16, and Kl 1 are seen to weave 



around the EEDD ring, making multiple contacts with residues in 
all four domains. We note that there is some redundancy here 
because two residues such as R13 and K16 can stiU cover all four 
domains consistent with the observation that //-GIIIA[K1 lA] can 
also block the channel [53]. The fact that two basic residues are 
sufficient to block the pore is also supported by several other )l- 
conotoxin blockers of Nay 1 channels, which have only two basic 
residues available at the binding interface (e.g., BuIIIB and SIIIA). 
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Figure 4. Bottom view of tlie complex structure demonstrating blocl<ing of the pore. The toxin residues R13, K1 1, and K16 mal<e contacts 
with the channel residues EEDD in all four domains. 
doi:1 0.1 371 /journal.pone.01 05300.g004 



Compared to potassium channels there are relatively fewer 
blockers of sodium channels, which is due to the larger pore size in 
the latter. A pore inserting Lys is sufficient to block a potassium 
channel whereas at least two basic residues are required to achieve 
the same in a sodium channel. The pore inserting Lys motif has 
been instrumental in functional studies of potassium channels 
using toxins peptides as probes. This has simplified interpretation 
of experimental results as well as construction of computational 
models of channel-toxin complexes. The situation in Nayl 
channels is much more complicated due to many possible 
configurations for coupling of 2-3 basic toxin residues with the 
EEDD residues in the pore. Also there are many high-affinity 
toxins that do not block Nayl channels. These features have 



certainly made interpretation of mutation experiments a more 
difficult task and sometimes resulted in conflicting proposals for 
the binding modes. Construction of accurate complex models 
using homology models of Nayl channels is expected to 
ameliorate this situation. Moreover such complex models will be 
very useful in designing analogues with enhanced affinity and 
selectivity properties, which may be required for development of 
toxin blockers of Nayl channels as therapeutic agents. 

Umbrella Sampling Simulations and Binding Free Energy 

We have previously shown in over a dozen case studies 
involving potassium channel toxins that the binding free energy 
can be determined near chemical accuracy from PMF calculations 
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Figure 5. Average rmsd of the /(-GIIIA backbone atoms for the residues 1 1-22 calculated for each umbrella window. The bull< value 
obtained from MD simulations of /i-GIIIA in a box of water is indicated by the dashed line. 
doi:1 0.1 371/journal.pone.01 05300.g005 



[58-63]. Thus calculation of the standard binding free energy of ;U- 
GIIIA will provide a complementary test for the accuracy of the 
proposed Navl.4— /i-GIIIA model. The PMF for the dissociation 
of /<-GIIIA from Nayl .4 is constructed from umbrella sampling 
MD simulation as described in Methods. Distortion of a ligand 
during umbrella ,samj)ling simulations is a concern in PMF 
calculations, especially for flexible peptides, which may lead to 
erroneous results if the distortion becomes permanent after the 
ligand [65] . We have checked against this possibility by calculating 
the average rmsd of /i-GIIIA in each umbrella sampling window. 
As shown in Figure 5, the toxin undergoes some distortion whUe it 
is puUed out of the binding pocket but the elastic energy associated 
with this distortion is recovered once the toxin is in the bulk region 
as indicated by the return of the rmsd to the bulk value. Lack of 
distortion is also verified by aligning the ^-GIIIA structures from 
the last umbrella window with the NMR structure. 

Two main questions in PMF calculations are how far the PMF 
should be extended and how long ('ach window should be run. 
The first is addressed by appearance of a flat region in the PMF 
which indicates that the ligand has reached the bulk region. The 
second question is rarely addressed in PMF calculations but it is 
equally important because without sufficient data, one is likely to 
obtain a wrong answer. W e address this issue by performing block 
data analysis of the PMF data. That is, we construct PMFs from 
2 ns blocks of data and sUde the blocks in 1 ns steps over the range 



of the available data. As shown in Figure 6, the PMFs initially 
drop monotonically and then fluctuate around a base line. During 
the first phase (1-4 ns), the PMFs drop because of the improved 
screening of the channel-toxin interactions as the system 
equilibrates. In the second phase (4-10 ns), fluctuations of the 
PMFs around a base line are of statistical nature and indicate that 
the system has been equilibrated. A common practice in PMF 
Ccdculations is to exclude an arbitrary amount of data for 
equilibration (e.g., 1 ns), and consider the rest of the data in 
production of the PMF. As seen from the convergence study in 
Figure 6, this could result in mixing of the equilibration and 
production data, which would lead to an overestimation of the 
binding afiinity. 

The fmal PMF for the Nayl .4-/i-GIIIA complex is indicated by 
a thick black line in Figure 6. We will comment on distinct features 
of the PMF when we discuss the evolution of the distances between 
the contact pairs below. Using Eq. 1 , we numerically integrate the 
PMF to find the binding constant and then determine the standard 
binding free energy using Eq. 2. The calculated value, —9.9 kcal/ 
mol, is in good agreement with the experimental value of — 10.5 
kcal/mol, which is determined from the most recent ICso 
measurement (19 uM) where state-of-the-art equipment was 
employed [74]. The level of agreement obtained for the standard 
binding free energy is similar to those obtained for potassium 
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Figure 6. Convergence study of the Nav1 .4-//-GIIIA PMF from block data analysis. To minimize fluctuations, a relatively large sampling size 
of 2 ns is used, which is slided in 1 ns steps over the 10 ns of data. The block-data PMFs drop monotonically in the first 4 ns as the system 
equilibrates and then fluctuate around a base line from 4-10 ns, signalling equilibration. The final PMF obtained from 4-10 ns is indicated by a thick 
black line. 
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channel toxins [58-64], and provides further support for the 
accuracy of the proposed model of the Navl.4-;U-GIIIA complex. 

The binding mode of the Navl.4— ^-GIIIA complex is 
dominated by charge interactions where the pairs are mostly at 
contact distances (Table 2). Analysis of how the contact distances 
change during dissociation provides complementary information 
on the relative strength of the various interactions as well as 
helping to explain specific features of the PMF. For this purpose, 
we have calculated the average pair distances in each umbrella 
window and plotted them as a function of the window position 
(Figure 7). The toxin residues are seen to break up in the following 
order and at the approximate positions: Kll at 30 A, K16 at 
34 A, R19 at 35 A, and finally R13 at 36 A. Because the weakest 
interactions will break up first and the strongest ones last, the 
persistence length of a pair during dissociation provides a good 
measure for its relative strength. This intuitive picture for the 
strength of the interactions is in good agreement with the mutation 
data in Table 2. For example, according to the set of data from 
[53], the affmity loss for mutations of Kl 1, K16, R19, and R13 to 
Ala are, respectively, 10, 81, 199, and 642. 

Apart from the glitch around 31-32 A, the final PMF in Fig 6 
follows a steadily rising trajectory until all the contact pairs are 



broken off at 36 A. The glitch appears to be associated with the 
fluctuations of the Kll and K16 pair distances (Figure 7). After 
z>36 A, the charged pairs interact via screened Coulomb 
interactions which corresponds to the shoulder region in the 
PMF. For z>40 A, the distances between the charged pairs are 
larger than 15 A. At those distances, the Coulomb interactions are 
mostly screened off, which correlates well with the PMF leveling 
off after z = 40 A and becoming flat for z>42 A (Figure 6). 

Conclusions 

In this work, we have constructed a homology model for the 
pore domain of the Nay 1.4 channel by aligning the DEKA 
residues in the selectivity filter with the corresponding EEEE 
residues in NayAb, whose crystal structure was determined 
recently. In order to validate the Nay 1.4 model, we have used 
the extensive functional data obtained from binding studies of /i- 
conotoxin GIIIA. An initial model for the Navl.4-^-GIIIA 
complex is created using HADDOCK, which is refined in MD 
simulations. The binding mode obtained is in broad agreement 
with the available mutagenesis data and shows that the toxin 
blocks the pore through multiple interactions of the R13, K16 and 
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Figure 7. Evolution of the distances between the interacting pairs as /t-GIIIA dissociates from Nav1.4. Average N-0 distances obtained 
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Kl 1 residues with the outer ring of EEDD residues in the channel 
vestibule. The standard binding free energy of //-GIIIA is 
determined from the PMF calculations and found to agree with 
the experimental value within chemical accuracy. Thus the 
proposed model of the Navl.4— /i-GIIIA complex has been well 
validated. Because there is a high degree of homology among the 
Nayl channels, the present Nay 1.4 model can be used as a 
template in constructing homology models for the pore domain of 
other Nayl channels. 

Our focus in this first study was the validation of the pore 
domain using the data on binding of //-GIIIA. For further studies 
of toxin binding to Nayl channels one needs to include the 
selectivity filter and the S5-P1 linker in the model. For example, to 
investigate the ionic strength dependence of toxin binding [54], a 
validated model of the selectivity filter is required. This can be 
achieved by studying the permeation and selectivity properties of 



Na""" ions, which we hope to tackle in a forthcoming paper. Our 
attempts to model the S,5-P1 linker in the turret of the Nay 1.4 
channel have not been successful due to lack of good templates. 
Because binding of /i-GIIIA does not appear to involve the S5-PI 
linker residues, a satisfactory binding mode could still be obtained 
without modeling this region. However, the S5-P1 linker residues 
are involved in binding of some other /i-conotoxins, and to 
understand the differences in their affinity and selectivity 
properties [74], it will be important to construct models of Nayl 
channels including the full turret region. The available mutagen- 
esis data could provide valuable guidance in this endeavor. Finally, 
there are ongoing efforts to harness the therapeutic potential of /H- 
conotoxins by designing analogues that selectively bind to a target 
Nayl channel [3-5]. Construction of accurate models of Nayl -/i- 
conotoxin complexes will be very useful in such efforts. 
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After the completion of this work, a paper on binding of ji- 
conotoxin PIIIA to the Nay 1.4 channel has appeared [75]. In this 
paper, two very different binding modes were predicted with 
almost equal binding free energies. This is wtry different from the 
unique binding mode found for ^-GIIIA in our study and needs to 
be investigated further. A snapshot comparing the Nay 1.4 models 
used in the two studies is given in Figure S I in File S 1 . 

Supporting Information 

File SI Figure SI, Snapshot comparing the pore domain of our 

Navl.4 model with that of Chen et al. [75]. 

(PDF) 
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